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ABSTRACT 

The  dot  diffusion  method  for  digital  halftoning  has  the 
advantage  of  parallelism  unlike  the  error  diffusion  method. 
However,  image  quality  offered  by  error  diffusion  is  still  re¬ 
garded  as  superior  to  other  known  methods.  In  a  recent  pa¬ 
per  we  showed  how  the  dot  diffusion  method  can  be  improved 
by  optimization  of  the  so-called  class  matrix.  In  this  paper 
we  first  review  the  dot  diffusion  algorithm  and  the  optimiza¬ 
tion  of  the  class  matrix.  A  method  for  inverse  halftoning 
of  dot  diffused  images  is  then  proposed.  The  method  uses 
wavelet  decomposition  to  eliminate  the  halftoning  noise  and 
does  not  make  use  of  the  knowledge  of  the  class  matrix. 


sion.  Since  dot  diffusion  also  offers  increased  parallelism,  it 
now  appears  to  be  an  attractive  alternative  to  error  diffusion. 

In  this  paper,  we  first  review  the  improved  dot  diffusion 
algorithm  of  [8],  and  then  address  the  inverse  halftoning 
problem.  Inverse  halftoning  has  a  wide  range  of  applications 
such  as  compression,  printed  image  processing,  scaling,  en¬ 
hancement,  etc.  In  these  applications,  operations  can  not  be 
done  on  the  halftone  image  directly,  and  inverse  halftoning 
is  mandatory.  A  simple  yet  efficient  algorithm  for  inverse 
halftoning  of  dot  diffused  images  is  proposed  and  compared 
to  other  methods. 

2  REVIEW  OF  DOT  DIFFUSION 


1  INTRODUCTION 

Digital  halftoning  is  the  rendition  of  continuous-tone  pictures 
on  displays  that  are  capable  of  producing  only  two  levels. 
There  are  many  good  methods  for  digital  halftoning:  ordered 
dither  [3],  error  diffusion  [4],  neural- net  based  methods  [2], 
and  more  recently  direct  binary  search  (DBS)  [10].  Ordered 
dithering  is  a  thresholding  of  the  continuous-tone  image  with 
a  spatially  periodic  screen  [3]  .  In  error  diffusion  [4],  the  error 
is  ‘diffused’  to  the  unprocessed  neighbor  points. 

Ordered  dithering  is  a  parallel  method,  requiring  only 
point  wise  comparisons.  But  the  resulting  halftones  suffer 
from  periodic  patterns.  On  the  other  hand  error  diffused 
halftones  do  not  suffer  from  periodicity  and  offer  blue  noise 
characteristic  [11]  which  is  found  to  be  desirable.  The  main 
drawback  is  that  error  diffusion  is  inherently  serial,  e.g.,  to 
get  the  halftoned  value  of  the  last  pixel,  all  of  the  remaining 
points  should  be  processed.  Also  there  occur  worm- like  pat¬ 
terns  in  near  mid-gray  regions  and  resulting  halftones  have 
ghosting  problem  [5].  Mitsa  and  Parker  have  optimized  or¬ 
dered  dither  matrix  [9]  for  large  size  like  256x256  to  get  the 
blue  noise  effect.  This  is  a  compromise  between  parallelism 
and  image  quality. 

The  dot  diffusion  method  for  halftoning,  introduced  by 
Knuth  [5],  is  an  attractive  method  which  attempts  to  retain 
the  good  features  of  error  diffusion  while  offering  substantial 
parallelism.  However,  surprisingly,  not  much  work  has  been 
done  on  optimization  of  the  so-called  class  matrix.  In  [8]  we 
showed  that  the  class  matrix  can  be  optimized  by  talking  into 
account  the  properties  of  human  visual  system  (HVS).  The 
resulting  halftones  are  of  the  same  quality  as  for  error  diffu- 
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The  dot  diffusion  method  for  halftoning  has  only  one  de¬ 
sign  parameter,  called  the  class  matrix  C.  It  determines 
the  order  in  which  the  pixels  are  halftoned.  Thus,  the  pixel 
positions  (ni,ri2)  of  an  image  are  divided  into  IJ  classes 
according  to  (t*i  mod  7,  ri2  mod  J)  where  I  and  J  are  con¬ 
stant  integers.  For  example,  Knuth  used  a  class  matrix  of 
size  I  =  J  =  8,  and  there  were  64  class  numbers  in  that  class 
matrix  [5].  Let  27(7*1,  712)  be  the  contone  image  with  pixel  val¬ 
ues  in  the  normalized  range  [0, 1].  Starting  from  class  A:  ==  1, 
we  process  the  pixels  for  increasing  values  of  k.  For  a  fixed 
fc,  we  take  all  pixel  locations  (711,7*2)  belonging  to  class  k 
and  define  the  halftone  pixels  to  be 


£h(7ii,n2) 


1  if  2(711,7*2)  >  0.5 
0  if  27(711,712)  <  0.5 


We  also  define  the  error  e(7*i,7*2)  =  27(7*1,7*2)  -  27^(7*1,712). 
We  then  look  at  the  eight  neighbors  of  (7*1,712)  and  replace 
each  contone  pixel  with  an  adjusted  version  for  those  neigh¬ 
bors  which  have  a  higher  class  number  (i.e.,  those  neighbors 
that  have  not  been  halftoned  yet).  To  be  specific,  neighbors 
of  x(i  j)  with  higher  class  numbers  are  replaced  with 


x(i,j)  -I-  2e(7ii,7i2)/u>  (for  orthogonal  neighbors)  (1  (a)) 

x(i,j)  +  e(7*i,  ri2)/iv  (for  diagonal  neighbors)  (1(b)) 

where  w  is  such  that  the  sum  of  errors  added  to  all  the 
neighbors  is  exactly  e(r*i,  7*2).  The  extra  factor  of  two  for  or¬ 
thogonal  neighbors  (i.e.,  vertically  and  horizontally  adjacent 
neighbors)  is  because  vertically  or  horizontally  oriented  error 
patterns  are  more  perceptible  than  diagonal  patterns. 

The  contone  pixels  27(7*1,7*2)  which  have  the  next  class 
number  k  +  1  are  then  similarly  processed.  The  pixel  values 
27(7*1,7*2)  are  of  course  not  the  original  contone  values  but 
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the  adjusted  values  according  to  earlier  diffusion  steps  (1). 
When  the  algorithm  terminates,  the  signal  (711,712)  is  the 
desired  halftone.  In  IJ  steps  the  algorithm  will  complete  the 
halftoning  process. 

Usually  an  image  is  enhanced  [5]  before  dot  diffusion 
is  applied.  For  this  the  continuous  image  pixels  C(i,j) 
are  replaced  by  C(i,j)  =  where  C(i,  j)  = 

£‘!i  C(i,j)/ 9.  If  a  =  0.9  then 

c\i,j)  =  8C(i,j)+C(i,j)-  ^  <?(«,«) ■ 

0<(u-*)2+(t>-j)2<3 

This  algorithm  is  completely  parallel  requiring  9  additions 
per  pixel,  and  no  multiplications. 

3  OPTIMIZATION  OF  CLASS  MATRIX 

Knuth  introduced  the  notion  of  barons  and  near  barons  in 
the  selection  of  his  class  matrix.  A  baron  has  only  low-class 
neighbors,  and  a  near-baron  has  one  high  class  neighbor. 
The  quantization  error  at  a  baron  cannot  be  distributed  to 
neighbors,  and  the  error  at  a  near  baron  can  be  distributed 
to  only  one  neighbor.  Knuth’s  idea  was  that  the  number 
of  barons  and  near  barons  should  therefore  be  minimized. 
He  exhibited  a  class  matrix  with  two  barons  and  two  near 
barons.  The  quality  of  the  resulting  halftones  are  still  infe¬ 
rior  to  error  diffusion  because  of  periodic  patterns  similar  to 
ordered  dither  methods  (see  Fig.  6).  In  our  experience,  the 
baron/near-baron  criterion  does  not  appear  to  be  the  right 
choice  for  optimization  as  explained  in  [8].  In  Sec.  3.1  we  in¬ 
troduce  a  different  optimization  criterion  based  on  the  HVS, 
and  show  that  the  image  quality  is  significantly  improved, 
though  the  class  matrix  does  not  minimize  barons. 

3.1  Objective  Function  Based  on  Blue  Noise 
It  has  been  observed  in  the  past  that  the  error  in  a  good 
halftone  should  have  the  blue  noise  property  [11].  This 
means  that  the  noise  energy  should  mostly  be  in  the  high 
frequency  region  where  it  is  known  to  be  less  perceptible.  We 
showed  in  [8]  how  to  incorporate  blue  noise  characteristics 
into  the  class  matrix  optimization: 

Let  Xh(ni,ri2)  denote  the  halftoned  version  of  a  constant 
gray  image  2(711,712)  =  g  where  0  <  g  <  l.Typically,  the 
dark  pixels  are  spatially  distributed  with  a  certain  aver¬ 
age  frequency  fg  called  the  principal  frequency,  which 
increases  with  gray  level  g.  Since  noise  energy  at  a  signifi¬ 
cantly  higher  spatial  frequency  than  fg  is  not  perceivable,  we 
can  optimize  a  halftoning  method  for  a  particular  gray  level 
g  by  forcing  the  noise  spectrum  to  be  concentrated  above  fg. 

Calculating  the  noise  spectrum.  In  order  to  implement  the 
optimization,  we  first  need  to  compute  the  noise  spectrum. 
The  halftone  pattern  2^(711, 712)  for  the  gray  level  2(711, 712)  = 
g  has  the  error  e(ni,ri2)  =  g  —  Xh (711,712),  which  is  an  N  x  N 
image.  As  explained  in  [8],  a  radially  averaged  power 
spectrum  Pr(kr)  for  this  error  is  calculated,  where  kr  is  a 
positive  integer  called  the  radial  frequency.  The  class  matrix 
in  the  dot  diffusion  method  should  be  optimized  such  that 
this  radial  spectrum  is  appropriately  shaped  for  a  well-chosen 
fixed  gray  level  g.  In  terms  of  the  radial  frequency  variable 
fcr,  the  principal  frequency  for  the  halftone  of  gray  level  g  is 
given  by 

fg  ~  kmaxy/g 


where  kmax  is  the  maximum  value  of  kr .  In  fact,  for  g  > 
0.5,  since  black  pixels  are  more  in  number,  the  halftone  is 
perceived  as  a  distribution  of  white  dots  and  we  have  to  take 

fg  —  km  ax  \/l  g  • 


Figure  1:  A  Typical  Figure  2:  The  Weight 
Desired  Radial  Spec-  Function  used  in  the  Op- 
tram  Characteristics  timization 

The  aim  of  the  optimization  is  to  shape  Pr(kr)  by  choice 
of  the  class  matrix  C  so  that  most  of  its  energy  is  moved 
to  the  region  kr  >  fg  (as  demonstrated  in  Fig.  1).  We 
therefore  define  the  cost  function 

rfg 

$(C,0)=  /  Pr(kr)w(kr)dkr 

Jo 

The  weight  function  was  chosen  to  be  w(kr)  =  (kr  —  fg)2 
for  0  <  kr  <  fg  and  zero  outside  [8].  In  the  optimization 
the  integral  was  replaced  with  a  discrete  sum.  The  choice 
of  the  class  matrix  that  minimizes  this  sum  was  performed 
using  the  pairwise  exchange  algorithm  [8],  which  was 
originally  proposed  in  digital  filter  literature  for  a  different 
application  [1].  The  gist  of  this  algorithm  is  as  follows: 

1)  Randomly  order  the  numbers  in  the  class  matrix. 

2)  List  all  possible  exchanges  of  class  numbers. 

3)  If  an  exchange  does  not  reduce  cost,  restore  the  pair  to 
original  positions  and  proceed  to  the  next  pair. 

4)  If  an  exchange  does  reduce  cost,  keep  it  and  restart  the 
enumeration  from  the  beginning. 

5)  Stop  searching  if  no  further  exchanges  reduce  cost. 

6)  Repeat  the  above  steps  a  fixed  number  of  times  and  keep 
the  best  of  the  class  matrix. 

Choice  of  gray  level.  Since  the  algorithm  can  be  applied 
only  to  a  given  gray  level,  the  gray  level  should  be  chosen 
wisely  to  get  good  halftones  for  other  gray  levels  also.  For 
most  natural  images,  the  best  gray  level  was  experimentally 
found  to  be  g  =  0.0625  as  explained  in  [8]. 

4  INVERSE  HALFTONING 

Inverse  halftoning  is  the  reconstruction  of  a  continuous  tone 
image  from  its  halftoned  version.  The  basic  aim  in  inverse 
halftoning  is  to  separate  the  halftoning  noise  from  the  orig¬ 
inal  image.  In  good  halftoning  algorithms,  the  noise  intro¬ 
duced  by  halftoning  is  concentrated  in  the  high  frequencies. 
Simple  low  pass  filtering  can  remove  the  high  frequency  noise 


but  it  also  removes  the  edge  information.  Thus  the  edge  in¬ 
formation  should  be  separated  from  the  halftoning  noise. 

Inverse  halftoning  using  wavelets  was  considered  in  [12] 
and  [6].  The  algorithm  in  [12]  is  tailored  for  error  diffusion, 
which  has  different  characteristics  than  dot  diffusion.  If  the 
method  in  [12]  is  used,  the  result  is  not  good.  This  can  be 
seen  from  Fig.  8  which  is  the  result  of  inverse  halftoning 
of  dot  diffused  Lena  by  using  the  method  in  [12].  The  im¬ 
age  suffers  from  periodic  patterns,  which  is  essentially  low 
frequency  noise. 

In  the  new  method,  the  specific  properties  of  the  dot  dif¬ 
fusion  algorithm  are  taken  into  account.  The  image  is  en¬ 
hanced  before  dot  diffusion,  hence  in  the  inverse  halftoning, 
the  dot  diffused  image  should  be  deenhanced  using  the  in¬ 
verse  filter  of  Fenh(zi,z2)  =  10  —  (z\  -fl  +  z^1)  (z2  +  1  + 22  x). 
Note  that  Fenh(e*Wl  >  0  for  all  0  <  w\ ,  W2  <  7r. 

We  use  the  wavelet  tree  built  from  the  analysis  block 
shown  in  Fig.  3.  An  image  C(x,y)  is  decomposed  into 
I/(z,y),  ff(x,y),  and  VXx,y)  using  the  undecimated  wavelet 
transform.  At  scale  2t+1,  (which  will  be  described  below), 
the  filtering  operations  are  as  follows  : 

L(wi,w2)  =  F(2twi)F{2tw2)C(wiyw2), 

H(wi,w2)  =  G(2lwi)F(2lw2)C(w1,w2), 

V(witw2)  =  F(2twi)G(2'w2)C(wi,w2)i 
where  G  and  F  are  derived  from  quadratic  spline  wavelets 
and  they  are  tabulated  with  the  synthesis  filters  in  Table  1, 
in  [7]  (F  is  H  in  the  latter  Table).  The  choice  of  filters 
given  in  [7]  detect  edges  at  different  scales  if  they  are  used 
in  the  wavelet  tree  shown  in  Fig.  4  with  scales  2°,21,22,23 
from  left  to  right.  For  example  ffj(x,y)  and  V%(x,y)  repre¬ 
sent  the  horizontal  edges,  and  vertical  edges  of  Li-i(x,y)  at 
scale  2*“*x  respectively,  and  Li(x,y)  is  the  low  pass  version 
of  Lt_i(x,  y). 


C(x,y) 


- >  L(x.y) 

- *  H(x.y) 

- >  V(x,y) 


Figure  3:  Wavelet  decomposition  of  an  image 
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Figure  4:  Wavelet  Tree  used  in  Inverse  Halftoning 


Let  us  denote  the  ith  level  low  pass  image,  vertical  edge 
image  and  horizontal  map  image  as  L*(x,  y),  Vi(x,y),  and 
Hi(x,y)  respectively.  The  4-level  wavelet  decomposition  is 


then  applied  to  the  deenhanced  halftone  image,  Lo(x,y). 
Then  for  each  pixel  location  (x,y),  the  following  is  done: 

1)  Apply  a  symmetric  FIR  Gaussian  filter,  /5(n,  m)  to 

Vi(x,y),  and  tfi(x,y).  (/ff(n,m)  =  ce  ^  for  -3  < 
n,  m  <  3,  and  c  is  chosen  such  that  the  DC  gain  of  the  fil¬ 
ter  is  unity).  The  first  level  edge  maps  contain  mostly  the 
halftoning  noise,  thus  low  pass  filtering  these  images  reduces 
the  blue  noise  without  harming  the  edges  too  much. 

2)  Let  E2z{x,y)  =  V2{x,y)Vz{x,y)  +  H2{x,y)Hz[x,y). 

if  £23(x,  y)  <  Ti  then  make  F2(x,  y)  =  0  and  H2(x,  y)  =  0. 

3)  Let  £34 (x,  y )  =  V3(x,y)Vi(x,y)  +  H3(x,y)H4(x,y). 

if  2?34(x,  y)  <  T2  then  make  V3(x,  y)  =  0  and  Hs(xi  y)  =  0. 

Steps  2  and  3  are  the  denoising  steps  in  the  algorithm.  In 
order  to  discriminate  the  edges  from  the  halftoning  noise,  we 
have  to  locate  the  edges.  For  this,  the  above  steps  perform 
a  cross  correlation  between  the  edges  at  different  scales.  If 
there  is  an  horizontal  edge  at  scale  i  at  (x,y)  then  /fi(x,y) 
and  f/i+i(x,y)  will  be  of  the  same  sign  [7].  The  same  is 
also  true  for  vertical  edges.  Combining  the  horizontal  and 
vertical  edge  correlations  gives  better  results  in  detecting  the 
diagonal  edges. 

4)  The  above  steps  have  modified  the  subband  signals  L», 
Hi  and  Vi  in  certain  ways.  We  now  use  the  inverse  filter 
bank  (synthesis  bank)  corresponding  to  Fig.  4,  and  obtain 
a  reconstructed  version  Lo(x,  y).  The  image  Lo(x,y)  is  the 
desired  inverse  halftone  image. 

5  EXPERIMENTAL  RESULTS 

The  512  x  512  continuous  tone  peppers  image  is  halftoned 
by  using  Knuth’s  class  matrix  (Fig  6),  and  by  the  optimized 
class  matrix  (Fig  7).  It  is  clear  that  the  new  method  is 
superior  to  unoptimized  dot  diffusion  method.  In  fact,  the 
new  method  offers  a  quality  comparable  to  FS  error  diffusion 
method  (Fig.  5).  Error  diffused  images  suffer  from  worm¬ 
like  patterns  which  are  not  in  the  original  image,  whereas 
dot  diffused  halftones  do  not  contain  these  artifacts.  Notice 
that  the  artificial  periodic  patterns  in  Fig.  6  are  absent  in 
Fig.  5  and  in  the  new  method  (Fig.  7). 

In  inverse  halftoning,  dot  diffusion  has  an  advantage,  even 
the  simple  unenhanced  image  is  a  quite  reasonable  inverse 
halftone  (psnr=26.62dB  for  Lena  image).  The  unenhanced 
image  is  further  processed  as  described  in  Sec.  4.  The  pa¬ 
rameters  used  in  the  method  are  found  experimentally.  The 
variance  of  the  Gaussian  filter,  a1  is  chosen  to  be  0.5  and 
the  thresholds  are  chosen  to  be  Ti  =  300  and  T2  =  20.  The 
results  are  shown  in  Fig.  10  (psnr=30.58dB)  and  in  Fig.  9 
(psnr=28.61dB).  2 

6  CONCLUSION 

Even  though  dot  diffusion  offers  more  parallelism  than  error 
diffusion,  it  has  not  received  much  attention.  This  is  partly 
because  the  noise  characteristics  of  error  diffusion  method 
are  generally  regarded  as  superior.  We  observed  that  by  op¬ 
timizing  the  class  matrix  for  blue  noise  at  a  fixed  gray  level, 
the  results  of  dot  diffusion  can  be  made  at  least  as  pleasing 

2  The  inverse  halftone  images  can  be  found  at 
htttp:/ / www.systems.caltech.edu  /  mese /halftone/ 


as  that  of  error  diffusion.  The  algorithm  terminates  in  at 
most  64  steps  for  8  x  8  class  matrix  compared  to  N 2  steps 
needed  for  error  diffusion  algorithm.  Moreover,  as  noticed 
in  [8],  the  algorithm  can  in  fact  be  terminated  in  about  50 
steps.  The  conclusion  is  that  Knuth’s  dot  diffusion  method 
with  a  carefully  optimized  class  matrix  is  very  promising; 
the  image  quality  is  comparable  to  error  diffusion,  and  the 
implementation  offers  more  parallelism  than  error  diffusion. 
We  also  developed  a  wavelet-based  inverse  halftoning  algo¬ 
rithm  which  works  very  well,  even  though  the  class  matrix 
information  is  not  used. 
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Figure  6:  Dot  Diffusion  with  Knuth’s  Class  Matrix 
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Figure  8:  Result  of  Inverse  Halftoning  using  previous 
method 


Figure  10:  Inverse  Halftoned  Peppers  Using  the  new 
method 


